
## Load packages ##
library(readxl)
library(haven)
library(tidyr)
library(dplyr)
library(ggplot2)
library(sp)
library(rworldmap)
library(countrycode)
library(RColorBrewer)

## Important information ##

#Tokelau is missing from rworldmap package
#Zanzibar does not have separate estimate

## Load data ##
data_hpacc <- read_excel("[paste path]/Data_map.xlsx")

## Generate ISO3 codes for matching ##
data_hpacc$country <- recode(data_hpacc$country, "England" = "United Kingdom", "STP" = "Sao Tome and Principe")
data_hpacc$country_iso <- countrycode(data_hpacc$country, "country.name", "iso3c")   

## Create dataframe for small islands nations ##
small_island <- data.frame(
  country = c("Barbados", "Bermuda", "Cabo Verde", "Comoros", 
              "Cook Islands", "Kiribati", "Malta",
              "Marshall Islands", "Nauru", "Niue", "Palau", "Saint Lucia", 
              "Samoa", "Seychelles", "Solomon Islands", "Sao Tome and Principe", 
              "Tuvalu", "Vanuatu", "Wallis and Futuna",
              "Tokelau", "Zanzibar"),
  lon = c(-59.5432, -64.7505, -23.8859, 43.3333, 
          -159.7771, 173.0000, 14.3754,
          171.1845, 166.9326, -169.8672, 134.5825, -60.9789, 
          -172.1046, 55.4919, 160.1562, 6.6131, 
          179.1945, 166.9326, -176.1946,
          -171.8182, 39.3167),
  lat = c(13.1939, 32.3078, 15.1201, -11.6455, 
          -21.2367, 1.8709, 35.9375,
          7.0644, -0.5477, -19.0544, 7.5149, 13.9094, 
          -13.7590, -4.6796, -9.6457, 0.1864, 
          -7.1095, -15.3767, -13.2825,
          -9.1668, -6.1333))
small_island <- small_island %>%
  left_join(data_hpacc, by=c("country")) 

## Define offsets manually ##
line_length <- 5  
manual_angles <- list(
  "Cook Islands" = 0, "Niue" = 320, "Samoa" = 0, "Wallis and Futuna" = 70,
  
  "Barbados" = 0, "Saint Lucia" = 40, "Bermuda" = 0,  
  
  "Cabo Verde" = 150, "Sao Tome and Principe" = 180, "Comoros" = 40, "Seychelles" = 0,
  
  "Vanuatu" = 30, "Tuvalu" = 0, "Solomon Islands" = 210, "Palau" = 0, "Nauru" = 180, 
  "Marshall Islands" = 30, "Kiribati" = 0,
  
  "Tokelau" = 70, "Zanzibar" = 20,
  
  "Malta" = 0
  )

## Convert angles to posiitons
deg_to_rad <- function(degrees) {
  return(degrees * pi / 180)
}
compute_offset_position <- function(country, lon, lat) {
  angle <- deg_to_rad(manual_angles[[country]])
  dx <- line_length * cos(angle)
  dy <- line_length * sin(angle)
  return(c(lon + dx, lat + dy))
}
adjusted_positions <- t(mapply(compute_offset_position, 
                               small_island$country, 
                               small_island$lon, 
                               small_island$lat))
small_island$lon_off <- adjusted_positions[, 1]
small_island$lat_off <- adjusted_positions[, 2]

## Adapt color palette ##
color_palette <- brewer.pal(8, "YlOrRd")  
get_color <- function(value) {
  bins <- seq(0, 80, length.out = 9)  # 8 categories (9 cut points)
  color_index <- findInterval(value, bins, rightmost.closed = TRUE)
  return(color_palette[color_index])
}
small_island$color <- sapply(small_island$glp1_elig_rel, get_color)
small_island$color[small_island$country == "Nauru"] <- "#8A001E"      # Manually color Nauru (beyond legend)
small_island$color[small_island$country == "Tokelau"] <- "#8A001E"      # Manually color Tokelau (beyond legend)

## Plot the world map ##
pdf("World map.pdf", width = 10, height = 7)
map_data <- joinCountryData2Map(data_hpacc, joinCode = "ISO3", nameJoinColumn = "country_iso", verbose = TRUE)
mapCountryData(map_data, 
               nameColumnToPlot = "glp1_elig_rel",
               mapTitle = "Proportion of adults (25-64 years) eligible for GLP-1 receptor agonists, %",
               colourPalette = color_palette, 
               missingCountryCol = "gray97",
               borderCol = "gray30",
               catMethod = seq(0, 80, length.out = 9))
segments(small_island$lon, small_island$lat, 
         small_island$lon_off, small_island$lat_off, 
         col = "black", lwd = 0.5)
points(small_island$lon_off, small_island$lat_off, lwd=0.1, pch = 21, bg = small_island$color, cex = 0.8)    
dev.off()






